#define NB_K2 10
#define _C2 0.05

template<class Kernel, class Resample>
class Particle<Kernel, Resample, Density::IBIS>{

public:
	Particle(Kernel *K,Resample *R,Density::IBIS *D,int M){
		_p=(*D).Get_p();
		mat Xt(M,_p);//p le nombre de parametre a obtenir de Density
		W=new double[M];
		X=Xt;
		_K=K;
		_R=R;
		_Z=0;
		_D=D;
		thres=_U*M;
		_M=M;
		this->Normalize();
		y=X.t();
		_n=(*D).Get_n();
		_K->Set_n(M);
	//	cout << "///" << _n << "///";
	}
	Particle<Kernel, Resample, Density::IBIS> (const Particle<Kernel, Resample, Density::IBIS>& X){
		cout << "Hello";
	}

	void Filter()
	{
		this->Init();
		for(int i=0;i<(_n);i++)
		{
		//	cout << i;
			cout << "step:" << i << "\n";
			this->Step(i);
		}
		
	}
	void Init()
	{
		//pior
		mat m(_p,1);
		mat foo(_p,_p);
		m.zeros();
		foo.fill(1);
		mat s=diagmat(foo);
		Distribution::Gaussian p(_p,m,s);
		//Distribution::Uniform p(_p,0,1);
		X=p.scrambled(_M);	
		//cout << X;	
		count_Skip=0;
	}
	
	void Step(int i)
	{
	
		/*if(count_Skip<_n-i-1)
		{
			double test=1;	
			while(test)
			{
				test=TestESS(i);
			}
		}*/
		double Ess=Ess_W();
		//cout << "Ess:" << Ess << "\n";
		cout << "Ess: " << Ess;
		if(Ess< thres)
		{
			cout << "Resample";
			mat v=growingvect(_M);
			double *w=new double[_M];
			for(int j=0;j<_M;j++){
				 w[j]=W[j];
			}
			(*_R)(&_M,w,v);
			delete[] w;
			//for(int i=0;i<_M;i++){cout << W[i] << "\n";}
			//cout << v << "\n";
			Arangemat(X,v);
			this->Normalize();//poids fixé a 1
			//_K->Set_s(X);
			_K->Set_bool(1);
			for(int t=0;t<NB_K2;t++)
			{
				_K->Move(&X,i);
				cout << "t " <<t;
			}
		}
		
		y=X.t();
				//dans le cas de tempering peut etre determiner
	//	cout << "before:" << y.col(1);
		_D->Weight(y,W,i);//le p designera comment on avance 
	//	cout << "after:" << y.col(1);
		X=y.t();
		cout << "mean : " << mean(X,0);
	}
	double Ess_W(double *w, int tf)//renormalisation???
	{
		double *foo= new double[_M];
		double sum=0;

		for(int i=0;i<_M;i++)
		{
		//	cout << w[i] << " ";
			sum+=w[i];
		}
		if(tf){
			_Z+=log(sum/_M);
		}
		for(int j=0;j<_M;j++)
		{
			foo[j]=(w[j]/sum);
		}

		sum=0;
		for(int i=0;i<_M;i++)
		{
			sum+=foo[i]*foo[i];
		}
		delete[] foo;

		return 1/sum;		
	}
	double Ess_W(void)//renormalisation???
	{
		double *foo= new double[_M];
		double sum=0;

		for(int i=0;i<_M;i++)
		{
			sum+=W[i];
		}
		_Z+=log(sum/_M);
		for(int j=0;j<_M;j++)
		{
			foo[j]=(W[j]/sum);
		}

		sum=0;
		for(int i=0;i<_M;i++)
		{
			sum+=foo[i]*foo[i];
		}
		delete[] foo;

		return 1/sum;		
	}
	
	mat Get_theta(void){
		return X;
	}

	void Normalize(void){
		for(int i=0;i<_M;i++){W[i]=1;}
	}
	
	mat   Get_W(void){
		double sum=0;
		mat C(_M,1);
		for(int i=0;i<_M;i++)
		{
			sum+=W[i];
		}	
		for(int i=0;i<_M;i++)
		{
			C(i,0)=(W[i]/sum)*_M;
		}
//		cout << C;
		return C;
	}
	int TestESS(int i)
	{
		double res;
		mat y=X.t();
		double *w=new double[_M];
		for(int j=0;j<_M;j++){
			 w[j]=W[j];
		}
		_D->Weight(y,w,i); 
		double ess=Ess_W(w,0);	
		delete[] w;
		if(ess<_C2*_M)
		{
			res=1;
			_D->Set_X(i);
			count_Skip++;
		//	cout << "Skip a data point";
		}else{
			res=0;
		}
		return res;
	}

	double Get_Z(void){
		cout << "ev: " << _Z;       
		return _Z;}
		////////////////////////////////////////////
	private:
	mat X;//value particle
	double *W;//weigths
	double Ess;
	int _M;
	double thres;
	double _Z;
	Kernel *_K;
	Resample *_R;
	Density::IBIS *_D;	
	mat y;
	int _n;
	int _p;
	int count_Skip;
};

